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Abstract 

This paper deals with the numerical study of a strongly anisotropic heat equation. The use of stan- 
dard schemes in this situation leads to poor results, due to the high anisotropy. Furthermore, the 
recently proposed Asymptotic-Preserving method [11 allows one to perform simulations regardless of 
the anisotropy strength but its application is limited to the case, where the anisotropy direction is given 
by a field with all field lines open. In this paper we introduce a new Asymptotic-Preserving method, 
which overcomes those limitations without any loss of precision or increase in the computational costs. 
The convergence of the method is shown to be independent of the anisotropy parameter < e < 1, and 
this for fixed coarse Cartesian grids and for variable anisotropy directions. The context of this work is 
magnetically confined fusion plasmas. 

Keywords: Anisotropic parabolic equation. Ill-conditioned problem. Singular Perturbation Model, 
Limit Model, Asymptotic Preserving scheme. Magnetic Island 



1. Introduction 

This work deals with the efficient numerical treatment of heat transport in a strongly anisotropic 
medium. We address in particular models of magnetised plasma with magnetic field perturbations such 
as those produced by tearing modes and magnetic islands. 

In classical transport theory of strongly magnetised plasmas, the ratio of the parallel (xy) to the 
perpendicular (x±)heat conductivity of a given species (electrons or ions) scales like (Q.c'TcY where Vtc 
is the cyclotron frequency (the rotation frequency around the field lines) and Tc the collision frequency. 
This product is several orders of magnitude (typically 10 to 12). 

Magnetic islands are non-ideal deformations of the primary magnetic field. In plasma confinement 
devices, they have a small magnetic component pointing outwards. However, due to the strong parallel 
conductivity, even a tiny outward components leads to a substantial heat loss in the island regions. 
Thus, magnetic islands are unwanted effects in actual applications. 

Theories of the formation of magnetic islands rely on various ingredients. In the regime where tearing 
modes (TM) are linearly unstable, magnetic islands are the result of TM evolution and saturation. 
When however TM are stable, magnetic islands can still occur through a mechanism of self-sustainment. 
In this regime, a key element of the island dynamics is the competition between the parallel and the 
perpendicular heat ffuxes, depending in particular on the ratio X||/x± (denoted in the sequel as l/s), 
which may ultimately determine whether the island grows or is suppressed. 
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The heat equation studied in this paper can be written as 

dtu - iV|| • (A|| V||ii) - • {A^V^u) = , (1) 

where Vy and Vx denote the gradient in the direction parahel (respectively perpendicular) to the 
magnetic field. This problem can become (depending on boundary conditions) ill-posed in the limit of 

Conventional numerical methods usually fail to give accurate results with reasonable computational 
resources for realistic physical parameters. Indeed, when l/e is of the order of 10^^ the problem becomes 
extremely anisotropic and standard discretizations lead to very badly conditioned linear systems with 
condition number proportional to l/{eh'^) {h being the spatial discretization step). It is therefore 
important to develop a numerical scheme that can address the problem and give accurate results 
independently of e. 

The original motivation of this work comes from the fusion plasma physics, but similar anisotropic 
problems are encountered in many other fields of application. One can mention for example image 
processing [161 [IS] ? transport modeling in fractured geological structures [2] or semiconductor modeling 

m- 

Numerical resolution of strongly anisotropic problems has been addressed by many authors. For 
example adapted coordinates are often used in the context of plasma simulation [ij [3l [9] . This approach 
can be however difficult to implement, especially when the magnetic field is variable in time. This is why 
it is preferable to choose a method which does not require mesh or coordinate adaptation, like in 
Another approach relies on numerical schemes specially developed for the anisotropic context. Finite 
difference schemes were investigated in [51 [171 [HI- High order finite element method was proposed 
in [7 . Multigrid methods [6 can sometimes be beneficial. These methods are usually efficient for a 
selected range of e but do not behave well in the limit e ^ 0. 

A different way to overcome this difficulty (adopted in this paper) is to apply the so called Asymp- 
totic Preserving scheme introduced first in [10 to deal with singularly perturbed kinetic models. The 
idea is to reformulate the initial problem into an equivalent form, which remains well-posed, even if the 
anisotropy strength is infinite. The reformulation that is studied in this paper was first applied to the 
anisotropic stationary diffusion equation in [5 and then to the nonlinear anisotropic heat equation in 
[T3l[TT]. This method is based on introduction of an auxiliary variable, which serves to eliminate from 
the equation the dominant part, i.e. the one multiplied by 1/s. The choice of the auxiliary variable 
presented in those papers allowed to solve the problem regardless of the anisotropy strength but im- 
posed serious limitations on the magnetic field. In particular, the case of magnetic islands cannot be 
treated by those schemes. In this paper we propose a new method which overcomes this limitation. 

The plan of the article is as follows. In Sec. [2] the mathematical problem is presented. Sec. [3] is 
devoted to the description of the numerical method. The numerical tests with known analytic solutions 
and the application to the problem of transport in a magnetic island are presented in Sec [4j 

2. Description of the mathematical problem 

We are interested in a resolution of an anisotropic, two or three dimensional heat problem defined 
on a domain Q. Let the anisotropy direction be given by a smooth and normalized vector field 6, \b\ = 1 
and let the computational domain be a bounded and sufficiently smooth two or three dimensional 
subset of with d = 2,3. The domain ft is equipped with a boundary F, which is decomposed 
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{PH) I 



accordingly to the boundary conditions into two parts: : Vd and Fjv = dQ.\VD with the Dirichlet and 
Neumann boundary condition imposed respectively. 

It is convenient to decompose vectors G M^, gradients V0, with (j){x) a scalar function, and 
divergences V • with v{x) a vector field , into a part parallel to the anisotropy direction b and a part 
perpendicular to it. These parts are defined as follows: 

v\\ := {v • b)b , v± := {Id — b0 b)v , such that v = v\\ -\- v± , 

V||0:= {b-V(l))b, V_L0:= {Id - b ®b)V(l) , such that V0 = Vy^ + V_l0 , 

V|| • t' := V • v\\ , V_L • t' := V • t'_L , such that V • v = V\\ • v -\- V_l • v , 

where (8) denotes the vector tensor product. 

The mathematical problem we are interested in reads: find the particle temperature u{t^ x), solution 
of the evolution equation 

f dtu- \Vw {A\\V\\u) - V ^- {A^V ^u) = , in [0,T]xl], 

\n\\ •(A||V||i^(t,-)) + ri^-(^±V^i^(t,-))=^iv(t,-), on [0,r] x T^v , 
u{t,')=gD{t,'), on [0,r]xrz^, 
u{{),-) = vP{-), in 

The diffusion coefficients A\\ and A±_ are bounded and of the same order of magnitude, satisfying 

< ^0 ^ ^11 (^) ^ ^1 : foi" almost all x G ^1, 
^oll'^lP < v^A^{x)v < Ai\\v\\^ , V'?; G M"^ and for almost ah x G 

with < Aq < Ai some positive constants. The anisotropy of the problem is characterized by a 
parameter 5, which can be very small and provoke substantial difliculties in the limit 5^0. 
Indeed, putting formally e = in (PH) leads to the following reduced problem 

-V||.(A||V||i/)=0, in [0,T]xl^, 

n||.(A||V||i.(t,.)) = 0, on [o,r]xr^,, 

u{t,')=gD{t,'), on [0,T]xrz^, 

1/(0,-) = u^{-), in n 

This problem may be ill-posed, depending on boundary conditions and the anisotropy field b. For 
example, when some field lines of b are closed in 17 the system would admit infinitely many solutions 
as any function constant along the closed lines of b (meaning V\\u = 0) and satisfies the boundary 
conditions solves the reduced problem. The same problem occurs when the field lines are open but do 
not pass through a boundary supplied with the Dirichlet conditions. This argument applies also to the 
case, where periodic conditions are imposed on a part of the boundary. This is the case in the numer- 
ical simulations related to the tokamak fusion plasma, where computational domain is topologically 
equivalent to a torus. Numerical discretization of the original (PH) problem in the limit £: ^ can 
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therefore lead to a very badly conditioned linear systems. In fact, the condition number is proportional 
to Ije. 

The goal of this work is to introduce a new numerical scheme which permits to solve the original 
singular perturbation problem (PH) independently of the value of e by the use of easy to implement 
numerical tools without a need of adaptation of any kind to the anisotropy strength/direction. The 
method presented herein generalizes the numerical scheme introduced in [11] . It removes the limitations 
of the cited numerical scheme while keeping all the advantages. It permits to solve the initial singular 
perturbation problem (PH) accurately on a simple Cartesian grid, which does not need to be specially 
tailored to suit the topology of the anisotropy field lines h and whose size does not depend on the 
anisotropy strength. 

The method is developed in the framework of Asymptotic-Preserving schemes. That is to say, 
the method is consistent with a so-called Limit problem as £ ^ and is stable independently of the 
small parameter e. The initial singular perturbation problem (PH) is reformulated in a suitable way. 
Introduction of an auxiliary variable allows to remove the terms proportional to 1/e from the equation. 
Resulting system is equivalent to the initial one and is well-posed in the limit £ ^ 0. The modification 
introduced here in allows to extend the applicability of the numerical scheme proposed in [TT] to the 
settings, where the field h may contain closed lines without any loss of accuracy. 

3. Numerical method 

S.l. Semi- discretization in space 



Let us write the variational formulation of the singular perturbation problem (PH) : find u{t^ •) G 
V := H^{Q) such that 



for almost every t G (0,T). As already discussed in the previous section, taking the formal limit of 
£ ^ leads to an ill-posed problem: 



with any function belonging to the vector space of functions constant in the anisotropy direction: 



being a solution. 

As proposed in [11 , a correct Limit problem can be established by seeking a solution in the subspace 
Q instead of V. In this case the leading order term (containing the parallel gradient) is eliminated from 
the equation and we are left with the following Limit problem: find u{t^-) G Q such that 




(2) 




a := {p G V / V||P = in 1^} 



(L) {dtu{t,-), v)v'y+ I 



A_\V ^u{t, ■) ■ V_lw dx 



L 
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for almost every t G (0,T). 

An Asymptotic Preserving scheme designed for the problem (P) should give accurate results and be 
stable independently of e, even in the limit s = 0. In particular, the condition number associated with 
corresponding linear system should not depend on e. That is to say, the correct Limit problem (L) has 
to be "hardcoded" in certain sense into the equations. Let us briefly recall the Asymptotic-Preserving 
scheme introduced in [11]. The method relies on the auxiliary variable q introduced by the relation 
€\/\\q = V||i^ in Q. This trick permits to get rid of the terms of order 0{l/e) from the variational 
formulation. The uniqueness of q is provided by setting g = on a part of the boundary Tin defined by 

Tin :={xeT / b{x) • n{x) < 0} , 

i.e. the part of the boundary, where the field lines enter the domain. The following reformulated 
problem, called in the sequel the Asymptotic-Preserving reformulation (AP-problem) is proposed: find 
■)^q{t, •)) G V X jCin^ solution of 



{AP) { 



(^,^)v*,v+ / {A±_V ±u) • V ±v dx ^ / A\\V\\q-V\\vdx - / gNvds = 0, 
Jn Jn Jfn 



A\\\/ \\u ■ \/ \\w dx — € / A\\\/\\q • \/\\w dx = 0^ \fw e Cin , 



yveV (3) 



where 



Cin := {q e L\Q) / Vyg G L\Q) and q\r,^ = 0}. 

The AP-problem is equivalent for fixed £ > to the original P-problem ([2|. Moreover, putting formally 
£ = in (AP) leads to a well-posed problem 



^^'^)^*'^+ / (A^V^ix) • V^vdx + / A\\V\\q-V\\vdx - f gNvds = 0, 



/ A\\V\\u ■ V\\w dx = 0^ Ww e Cin , 



which is equivalent to the correct Limit problem (L). In this case, the auxiliary variable q acts as a 
Lagrange multiplier forcing u to be constant along b. 

The drawback of this method is the choice of the space for the auxiliary variable. Imposing ^Ir^^ = 
provides uniqueness of a solution but limits the application of the scheme to the case where all field 
lines are open. Indeed, fixing a value of q on the infiow boundary does not provide uniqueness of q 
on field lines which does not intersect with the infiow boundary {i.e. on closed field lines). In order 
to overcome this restriction we propose a new approach based on penalty stabilization rather than on 
fixing the value of q on one of the boundaries. The modification of the second equation of the AP scheme 
^ consists of an introduction of a penalty term — a mass matrix qw multiplied by a stabilization 
constant. A suitable choice of this constant permits to conserve an accuracy of the scheme. 
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The new method (APS-scheme) reads: find {u{t^ ')^Q{t^ •)) G V x £, solution of 



^ ^^'^)v%v+ / {A^V±u) -V^vdx ^ I A\\V\\q-V\\vdx - [ gNvds = 0, 
Jn Jn Jfn 



(APS) < 



where 



A\\\/ \\u ■ \/ \\w dx — £ / A\\\/ \\q ■ \/ \\w dx — a / qw = 0^ \/w ^ C ^ 
n Jn Jn 



yveV 



(4) 



C:={qe L^{Q) / V\\qe L^{9)} 



and a is a positive stabilization constant. We postpone a theoretical justification of the proposed 
method to the forthcoming paper. A similar method with stabilization by a diffusion matrix instead of 
a mass matrix was recently proposed in the context of a a posteriori error indicator and mesh adaptation 
for strongly anisotropic elliptic equations in [14 . 

Let us now choose a polygonalization of the domain ^ with polygons of the diameter approximately 
equal to h and introduce the finite element spaces V/i C V and Ch C The finite element discretization 
of Q writes then: find (uh^qh) ^ Vh x Ch such that 



{APS)h I 



I ^^Vhdx^ I {A^V^Uh)-V^Vhdx^ [ A\\V\\qh ■ V\\Vhdx - [ gNVhds = 0, 
Jn Jn Jn Jfn 

/ A\\\/\\Uh 'S/\\Whdx - £ / A\\\/\\qh -S/wWhdx - h^^'^ / qw = 0, \/w e Ch 
< Jn Jn Jn 



yvh e Vh 



(5) 

Remark that in order to ensure convergence rate in L^-norm we have put a = where k is the 

order of finite element method. 

3.2. Semi- discretization in time 

One should be extremely careful when discretizing in time the ([5| scheme as not all numerical 
schemes conserve the AP property. In fact a chosen method should be L-stable. This is the case for 
a standard first order implicit Euler scheme. If however a higher order in time method is desired then 
certain simple schemes are excluded. For example a Crank-Nicolson discretization is only A-stable and 
not L-stable. The obtained system would not be asymptotic preserving giving reliable results only 
under certain assumptions. Namely, £ should be close to one, time step should be of the order of £ or 
the initial value should already be a solution to the stationary equation, i.e. the parallel gradient of 

should be proportional to £. This is why we choose in this work to present both a first order implicit 
Euler scheme and a second order, L-stable Runge-Kutta method. 



3.2.1. Implicit Euler scheme 
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Let us introduce the standard bilinear forms 

Jn 

«||(©,X):= / A||V||e. Vyx^ix, a^(e,x):= / A^V • V dx , 
Jn Jn 

and write the first order, implicit Euler method in more compact notation: Find {u^^^ ^Q^^^) ^ ^h^^h^ 
solution of 



{Ea. 



PS) 



{ul+\vh)+T(a^{ul+\vh)+a^^{q]:+\vh)-J^^gN{t^+')vas) = 
a^^(ul+\wh)-ea^^{ql+\wh) - h''+\ql+\wh) =0, 



3.2.2. L-stable Runge-Kutta method 

Any s-stage Runge-Kutta method applied to following problem 



du 
'di 



is conveniently defined using a Butcher diagram: 





ail • 


• CLis 


Cs 




^ss 




bi • 





The method reads: for given u'^, an approximation of u{tn), the u'^^^ is a linear combination the 
method 



,n+l 



where Ui are solutions to the following problems: 



Ui = u'^ -\-T^^aij{Luj + f{t + Cjt)), . 

Moreover, if bj = agj for j = 1, . . . , 5 than u^^^ is equal to the last stage of the method: u^^^ = Ug. 

In order to obtain a second order accurate in time scheme, we choose to implement a two stage 
Diagonally Implicit Runge-Kutta (DIRK) second order scheme. The scheme is developed according to 
the following Butcher's diagram: 



A 


A 





1 


1-A 


A 




1-A 


A 


= 1 


1 

V2- 
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The second order AP-scheme writes: find {u^'^^ , q]^'^^) G Vh x solution of 



.n+1 _ n+1 



^n+1 _ ^n+1 



with u^\^ and 1^2^^ denote the solutions of the first and the second stage of the Runge-Kutta method. 



4. Numerical results 

We are now ready to test the proposed schemes. Let us perform in the beginning numerical ex- 
periments for b having all field lines open. This setting allows us to compare implicit Euler-APS and 
DIRK- APS methods with their non penalty stabilized versions (Euler-AP and DIRK-AP) proposed in 
pT] . The Euler-AP and DIRK-AP schemes are obtained by discretization of the AP-problem and differ 
from Euler-APS and DIRK- APS in two details. Stabilization terms are absent and the vector space 
is replaced by jCin,h- The methods are also compared to a standard implicit Euler discretization of the 
initial singular perturbation problem (PH), given by 

{P)„r {ul+\vh) + T(^a^{ul+\vh) + \a\\r,i{ul,ul+\vh)- QNit^'+^Hds^ ={utv„). 

Finally, the DIRK- APS scheme is applied to the case of magnetic islands, where some field lines of 
b are closed. But let us first introduce a finite element space that we intend to use in all numerical 
experiments. 

4-1- Discretization 

Let us consider a 2D square computational domain Q = [0,1] x [0,1]. We choose to perform all 
simulations on structured meshes. Let us introduce the Cartesian, homogeneous grid 

Xi=i/N:, , 0<i<N^, yj=j/Ny , 0<j<Ny, 

with Nx and Ny being positive even constants, corresponding to the number of discretization intervals 
in the x- resp. ^-direction. The corresponding mesh-sizes are denoted by hx > resp. hy > 0. We 
choose a standard Q2 finite element method (Q2-FEM), based on the following quadratic base functions 
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' (x-Xi-2)ix-Xi-i) 
/ {Xi + 2-x){Xi + i-x) 





X G [Xi,Xi^2] 

else 



(y-yj-2)(y-yj-i) 
2hl 

iyj+2-y)iyj+i-y) 
2hl 





for even z,j and 

f {XiJ^i-x){x-Xi-i) 



hi 







X e [xi-i,Xi^i\ 

else 



-y)iy-yj-i) 



for odd j, we define the space 



W/, := {vh = Vij 6>^^ {x) Oy^ {y)} . 



1,0 



The spaces Vh and Ch are then defined by 

Vh = jC^h = Wh, Cin^h = {qh ^ , such that (7/,|r,, 



y e [yj-2,yj], 
y e [yj,yj+2], 
else 



y ^ [%-i.%+i] 
else 



0}. 



The matrix elements are computed using the 2D Gauss quadrature formula, with 3 points in the x and 
y direction: 



where Co = = 0, C±i = V±i = ^yl' ~ ^/^ '^±i =5/9, which is exact for polynomials of 
degree 5. Linear systems obtained for all methods in these numerical experiments are solved using a 
LU decomposition, implemented by the MUMPS library. 

4-2. Numerical tests 

4-2.1. Known analytical solution 

Let the computational domain ft be supplied with the boundaries Tn = {(0,x) U G [0, 1]} 

and Td = {(x, 0) U (x, l)\x G [0, 1]} with the boundary conditions Qn — go = 0. 

Let us now construct a numerical test case with a known solution. Finding an analytical solution 
for an arbitrary 6-field is extremely difficulty. We prefer rather to act as in the previous papers [H [TTj , 
where we started from selecting a suitable limit solution uq and than defined anisotropy direction as 
its isolines. Let 

uo = sin [iry + a{y'^ — y) cos(7rx)) e~^, 

where a is an amplitude of variations of h. For a = 0, the limit solution is constant in the X-direction 
and represents a solution for h = (1,0)^. For a 7^ the anisotropy direction is determined by the 
following implication 



bx 



dx 



dy 



0, 
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h 


L^-error e = 1 


P 


Eap 


Eaps 


RKap 


RKaps 


0.1 


5.6 X 10-^ 


5.6 X 10-^ 


5.6 X 10-^ 


5.6 X 10-^ 


5.6 X 10-^ 


0.05 


7.1 X 10-"^ 


7.1 X 10-"^ 


7.1 X 10-"^ 


7.1 X 10-"^ 


7.1 X 10-"^ 


0.025 


8.9 X 10-^ 


8.9 X 10-^ 


8.9 X 10-^ 


8.9 X 10-^ 


8.9 X 10-^ 


0.0125 


1.11 X 10-^ 


1.11 X 10-^ 


1.11 X 10-^ 


1.11 X 10-^ 


1.11 X 10-^ 


0.00625 


1.39 X 10-^ 


1.39 X 10-^ 


1.39 X 10-^ 


1.39 X 10-^ 


1.39 X 10-^ 


0.003125 


1.74 X 10-'^ 


1.74 X 10-'^ 


1.74 X 10-'^ 


1.74 X 10-'^ 


1.74 X 10-'^ 


h 


L^-error e = 10"^° 


P 


Eap 


Eaps 


RKap 


RKaps 


0.1 


6.9 X 10-1 


1.62 X 10-^ 


1.66 X 10-^ 


1.62 X 10-^ 


1.66 X 10-^ 


0.05 


6.9 X 10-1 


2.20 X 10-"^ 


2.37 X 10-"^ 


2.20 X 10-"^ 


2.37 X 10-"^ 


0.025 


6.9 X 10-1 


2.77 X 10-^ 


2.93 X 10-^ 


2.77 X 10-^ 


2.93 X 10-^ 


0.0125 


6.9 X 10-1 


3.43 X 10-^ 


3.58 X 10-^ 


3.43 X 10-^ 


3.58 X 10-^ 


0.00625 


6.9 X 10-1 


4.2 X 10-'^ 


4.4 X lO-'^ 


4.2 X 10-'^ 


4.4 X lO-'^ 


0.003125 


6.9 X 10-1 


5.3 X 10-^ 


5.5 X 10-^ 


5.3 X 10-^ 


5.5 X 10-^ 



Table 1: The absolute error of u in the L^-norm for different mesh sizes and e = 1 or e = 10-^^, using 
the singular perturbation scheme (P) and the two proposed AP-schemes for a time step of r = 10-^ 

and at instant t = 10-^. 



and therefore h can be defined for example as 

^ _ B — ( ^'^^^ ~ COs(7rx) + TT \ 

~ |5| ' ~ \ 7ra(?/^ - ?/) sin(7rx) J 

We set a = 1 in all our simulations so that the direction of the anisotropy is variable in the computational 
domain. Note that we have 5 7^ in the computational domain. 

Finally, a perturbation proportional to e is added to and an a function u is obtained. In this 
way we ensure that u converges, as £ ^ 0, to the limit solution uq. For example 

u = sin (ny + a{y'^ — y) cos(7rx)) + e cos {2tix) sin (ny) . (6) 

We set the initial condition to be equal to u{t = 0), with u defined by ^ and add a suitable force 
term to the right hand side of the numerical schemes so that u is an exact solution to the problem. In 
this setting we expect all Asymptotic-Preserving methods to converge in the optimal rate, independently 
on e. 

In order to validate our method and confirm our expectations we verify first the convergence in 
regard of the space discretization. We choose a time step small enough so that the space discretization 
error is much bigger than the time discretization error. Then, we perform numerical simulations for 
100 time steps for different mesh sizes for all Asymptotic-Preserving schemes and the implicit Euler 
discretization of the (P) problem. Numerical errors are given in Table [l] and Figures [l] and ^ The 
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third order space convergence in the L2-norm for large values of e is achieved (as expected) by all 
methods. Stabilization procedure does not alter the accuracy for weak anisotropy. For small values of 
e only the Asymptotic Preserving schemes give good numerical solutions. In this case the stabilization 
term decreases slightly the precision of the results (by 2.9 — 4.6% compared to non stabilized schemes), 
keeping however the order of convergence. 



1e-05 



1e-06 



1e.07 



0.001 




1 

0.1 
0.01 
0.001 
0.0001 
1e-05 
le.06 
1e^07 



le.08 



0.001 



' (AP-E 
(APS-E) 
(AP-RK) 

(APS-RK) 



0.01 0.1 

-20 



(b) £ = 10- 



Figure 1: Relative L^-errors between the exact solution and the computed solution for the standard 
scheme (P), the Euler-AP method {Eap) the stabilized Euler-APS method {Eaps), the DIRK-AP 
scheme (RKap) and the stabilized DIRK- APS scheme (RKaps) as a function of /i, for £ = 1 resp. 

e = 10~^^ and the time step r = 10~^. Observe that for e = 1 all schemes give the same precision, for 
e = 10~^^ the standard scheme does not work while all AP schemes give comparable accuracy. 



Next, the time convergence of the methods is tested. In this case the mesh size is chosen in such 
a way that the space discretization does not influence the numerical precision (at least for large time 
steps). We perform simulations for different time steps and for a fixed final time {t = 0.1). The results 
are presented in Table [2] and Figure |3] The time discretization convergence order is confirmed. The 
{RKap) and {RKaps) schemes are, as expected, of second order in time as long as the error due 
to the space discretization is smaller than the error induced by the time discretization. The {Eap) 
and (Eaps) schemes are of first order for all values of the anisotropy parameter and the standard 
(P)-scheme works well and is of first order where the anistropy is weak {e close to one). Note that the 
stabilization procedure does not influence the accuracy of the solution if the time discretization error 
is bigger than the space discretization error. One can also observe that reasonably accurate results are 
obtained even with relatively large time steps, especially for the Runge-Kutta schemes. This property 
would allow orders of magnitude gains in integration time for a given physics application with respect 
to non-AP schemes. 

Numerical test have confirmed that all asymptotic preserving schemes are accurate independently 
of e. The penalty stabilization procedure can introduce a very small amount of additional error in some 
configurations, but the order of convergence is conserved. In the following experiments, the stabilized 
schemes are tested on the case of b containing closed field lines. 
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T 


L^-error s = 1 


P 


Eap 


Eaps 


RKap 


RKaps 


0.1 


1.32 X 10-^ 


1.32 X 10-^ 


1.32 X 10-^ 


8.4 X 10-5 


8.4 X 10-5 


0.05 


7.2 X 10-^ 


7.2 X 10-^ 


7.2 X 10-^ 


2.43 X 10-5 


2.43 X 10-5 


0.025 


3.55 X 10-^ 


3.55 X 10-^ 


3.55 X 10-^ 


6.2 X 10-^ 


6.2 X 10-^ 


0.0125 


1.75 X 10-"^ 


1.75 X 10-"^ 


1.75 X 10-"^ 


1.59 X 10-^ 


1.59 X 10-^ 


0.00625 


8.8 X 10-^ 


8.8 X 10-^ 


8.8 X 10-5 


4.8 X 10-^ 


4.8 X 10-'^ 


0.003125 


4.4 X 10-^ 


4.4 X 10-^ 


4.4 X 10-^ 


2.82 X 10-5 


2.82 X 10-5 


0.0015625 


2.20 X 10-^ 


2.20 X 10-^ 


2.20 X 10-5 


2.64 X 10-^ 


2.64 X 10-^ 


T 


L^-error s = 10"^° 


P 


Eap 


Eaps 


RKap 


RKaps 


0.1 


2.28 X 10-1 


1.14 X 10-^ 


1.14 X 10-3 


7.4 X 10-5 


7 A X 10-5 


0.05 


2.53 X 10-1 


6.2 X 10-"^ 


6.2 X 10-"^ 


2.10 X 10-5 


2.10 X 10-5 


0.025 


2.53 X 10-1 


3.07 X 10-^ 


3.07 X 10-"^ 


5.3 X 10-^ 


5.3 X 10-^ 


0.0125 


2.50 X 10-1 


1.51 X 10-"^ 


1.51 X 10-"^ 


1.33 X 10-^ 


1.33 X 10-^ 


0.00625 


2.51 X 10-1 


7.6 X 10-5 


7.6 X 10-5 


3.44 X 10-^ 


3.44 X 10-^ 


0.003125 


2.53 X 10-1 


3.81 X 10-5 


3.81 X 10-5 


1.18 X 10-^ 


1.18 X 10-^ 


0.0015625 


2.53 X 10-1 


1.90 X 10-5 


1.90 X 10-5 


8.3 X 10-^ 


8.3 X 10-^ 



Table 2: The absolute error of u in the L^-norm for different time step using the singular perturbation 
scheme (P) and two proposed AP-schemes for mesh size 200 x 200 at time t = 0.1. 



12 



1 

0.1 
0.01 
0.001 
0.0001 
le-05 
^e-06 
le.07 



- (P) 


1 


A (AP-E) 




\ (APS-E) — •— 


0.1 


\ (AP-RK) 


XAP&RK) — y 






0.01 




0.001 




0.0001 




le.05 




1e-06 




le.07 



1e-15 le.10 

(a) h = 0.1 
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Figure 2: Relative -errors between the exact solution and the computed solution for the standard 
scheme (P), the Euler-AP method {Eap) the stabilized Euler-APS method {Eaps), the DIRK-AP 
scheme {RKap) and the stabilized DIRK-APS scheme {RKaps) as a function of for h = 0.1 resp. 
h = 0.00625 and the time step r = 10~^. Observe that the standard scheme is accurate only for small 
values of anisotropy (large e) while all AP schemes give comparable accuracy in the whole range of e. 
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Figure 3: Relative -errors between the exact solution and the computed solution with the 
standard scheme (P), the Euler-AP method (Eap) the stabilized Euler-APS method (Eaps)^ the 
DIRK-AP scheme {RKap) and the stabilized DIRK-APS scheme {RKaps) as a function of r, and for 
e = 1 resp. e = 10~^^ and a mesh with 200 x 200 points. 



4.3. Temperature balance in the presence of magnetic islands 

In this section we perform a numerical experiment related to the tokamak plasma. This numerical 
test case fully demonstrates the novelty of the proposed stabilized AP scheme as it can applied in more 
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general settings as the previous ones [5l [TT]. The results of this section show that it is capable of 
simulating the heat transfer in the presence of the closed field lines in the computational domain. 

We consider a square computational domain Q = [—0.5,0.5] x [—0.5,0.5] and a field b with a 
perturbation consisting of a region with closed lines. This so called magnetic island is initially localized 
in the center of the domain and moving with the velocity uj. The field is given by 



where A is a perturbation parameter related to the island's width w = AA^^'^/n. This is the largest 
distance between the two branches of the separatrix, the line that divides the domain into regions of 
open and closed field lines. The two branches meet at the X-point, the saddle point of the vector 
potential. The island center, an extremum of the vector potential, is referred to as the 0-point. If 
A = the obtained field is aligned with the Y axis and points upwards (downwards) for x > (x < 0). 
For A>0 the magnetic island consisting of closed field lines appears in the region around x = 0. 

Referring to the theory of magnetic islands, this magnetic field geometry approximates a saturated 
tearing mode in the so called constant-?/^ regime (whence the parameter A is a constant). The frequency 
describes the rotation of an island as observed in experiments and numerical simulations. 

This frequency is typically of the order of the diamagnetic frequency, which is smaller than the 
Alfven frequency (the propagation rate of an Alfven wave), but larger than the transport rate across a 
typical island. It is thus interesting to study heat transport in an island that rotates sufficiently fast. 

For a static island or one rotating sufficiently slowly, we expect the fast transport along the field 
lines to fiatten the temperature profile in the island region. 

We choose to impose periodic boundary conditions on {{x^y) G dQ \y = —0.5} U {{x^y) G \y = 
0.5} so the domain is topologically equivalent to the surface of a torus. We supply the computational 
domain with two sets of conditions on remaining boundaries: 

1. Dirichlet boundary conditions Td = {{x^y) G dft \x = —0.5} U {{x^y) G \x = 0.5} with 



where temperature is exchanged with the exterior only by diTij boundary. 

2. Neumann and Dirichlet boundary conditions: Tn = {{x^y) G \x = —0.5} and Tjj = {{x^y) G 
dn \x 0.5} with 



which corresponds to the constant heating of the left side of the computational domain. 

In the case of the boundary conditions of the first type we expect that the presence of the island should 
increase the gradient of the temperature outside the island region and keep the temperature constant 
inside the island in such way that the total energy of the system remains unchanged. If the boundary 
conditions are of the second type, i.e. in the heating case, the gradient of the temperature in the 
non-island region should remain constant leading thus to a loss of the total energy of the system. 

We perform simulations with fixed e = 10~^^ and A = 0.01, giving an island of a width w = 
0.4:/ TT ~ 13%. The magnetic field lines are shown on Figure [4j The initial conditions for both boundary 





for X = —0.5 
for X = 0.5 



gD{tr) = 0, 



(7) 
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Figure 4: Magnetic island for A = 0.01 



types are the same and correspond to the stationary solution with no island present. That is to say, 
u^{x^y) = —X + ^. We perform our simulations on a fixed grid of 200 x 200 points for 100 time 
steps r = 2.5 x 10~^ until the final time 0.25 is reached. We compare results for a stationary island 
{uo = 0) and a moving island {uo = 10) with no island case. We are interested in the temperature profile 
along the X axis as well as in the total energy of the system, i.e. integral of the temperature in the 
computational domain. 

4.3.1. Dirichlet boundary condition on the left edge ^ 

Numerical results confirm our expectations. The total energy remains constant in the system. 
Integral of the temperature in the computational domain equals to 1/2 in all three cases: for a stationary 
island, moving island and non perturbed system. The temperature is constant in the island region 
leading to a stronger gradient outside the perturbation. Temperature profiles along the X axis are 
shown on Figures [5] and [6] In the case of a moving island the width of a constant temperature region 
in the temperature profile along the X axis is oscillating in time as the island moves in the domain. It 
is interesting to note that even at the time when there are no closed field lines across the X axis the 
fiat region can be observed. In fact, the x component of temperature gradient is negative but close to 
zero in this region. Remark also that the temperature gradient increases, approaching zero, near the 
island at its largest width, i.e. across the 0-point. One refers to this effect as "profile fiattening". 

4.3.2. Neumann boundary condition on the left edge ^ 

In the case of Neumann boundary condition imposed on the left boundary of the domain the results 
again are consistent with our expectation. The presence of the island reduces the total energy. Integral 
of the temperature drops from 0.5 to 0.44 and the maximal temperature from 1 to 0.89 for both 
stationary and moving island. Temperature profiles along the X axis are shown on Figures [7| and 
|8] The results are very similar to those in the previous test case. The difference is in the maximal 
temperature and in the temperature gradient in the X direction, which is now close to one in the 
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Figure 5: Temperature profiles (top row) and temperature gradients (bottom row) along the X axis 
for non perturbed field {A = 0) on the left and a stationary island {A = 0.01) present in the center of 

the domain on the right for the Dirichlet BC. 



regions far from the island. This is of course what is expected and consistent with the heating imposed 
on the left boundary by the means of Neumann boundary condition. 

4.3.3. Fast rotating magnetic island with Dirichlet boundary condition on the left edge ^ 

In the last experiment we study the effect of the islands rotation speed on the temperature profile. 
We expect to see a different temperature profile when the rotation is faster than the transport rate in 
the parallel direction. To achieve this numerically we decrease the island width such that A = 0.000625, 
augment the mesh size to 500 x 500, decrease the size of the computational domain Q = [—0.125, 0.125] x 
[—0.5,0.5] and put e = 10~^. We also increase the time resolution and put r = 0.25 x 10~^ and vary 
the rotation speed from 10^ to 10^. Then we compare temperature profiles after 10000 time steps. 

For the smallest rotation velocity {uj = 10^) no deviation from the stationary case is observed. For 
other velocities the temperature profile across the islands center is only slightly affected by the rotation 
velocity. The most interesting effect of the velocity is visible on the profile taken away from the island. 
In the stationary case the profile is a straight line with a constant gradient. If the rotation velocity is 
big enough, the profile begins to slightly flatten near the center (X = 0) for uj = 10^ and finally looks 
exactly the same as for the islands center (a; = 10^ and cu = 10^). In fact, for a sufficiently big rotation 
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Figure 6: Temperature profiles along the X axis for a moving island {A = 0.01, uj = 10) in the first 
row, X component of temperature gradient in the second row and a corresponding anisotropy field in 
the last row for different time steps for the Dirichlet BC. 



speed the temperature profile becomes homogeneous, i.e. independent of Y. A zoom of temperature 
profiles at the islands center {Y = Yc) and at the most distant from the islands center Y {Y = 1^ + 0.5 
if Yc < and Y = Yc — 0.5 otherwise) is presented on Figure |9] 
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Figure 7: Temperature profiles (top row) and temperature gradients (bottom row) along the X axis 
for non perturbed field {A = 0) on the left and a stationary island {A = 0.01) present in the center of 

the domain on the right for the Neumann BC. 

5. Conclusion 

The here presented Asymptotic-Preserving scheme proves to be an efficient, general and easy to 
implement numerical method for solving strongly anisotropic parabolic problems. As numerical exper- 
iments show, the second order two stage Runge-Kutta scheme is of particular interest as it produces 
accurate results for relatively large time steps allowing to substantially reduce the computational time. 
The herein introduced stabilization term is an important improvement to the previously proposed 
Asymptotic Preserving schemes. The not so rare in real applications anisotropy topologies can be 
successfully addressed with this new method. Finally, the case of magnetic islands studied in the last 
section brings new and interesting results showing the homogenization of the temperature profile for 
fast rotating islands. 
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Figure 9: Comparison between temperature profiles along the X axis for the stationary and rotating 
islands. The profile across the island center (Y = Yc) on the left and away from the center 

{Y = Yc± 0.5) on the right. 
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